﻿



#include <stdio.h>
#include <string.h>
#include <assert.h>
#include "XMen_sm2.h"
#include "XMen_asn1.h"


/**
 * @brief SM2 prime field modulus p
 */
const SM2_BN SM2_P = {
	0xffffffffu, 0xffffffffu, 0x00000000u, 0xffffffffu,
	0xffffffffu, 0xffffffffu, 0xffffffffu, 0xfffffffeu,
};


/**
 * @brief SM2 curve parameter b
 */
const SM2_BN SM2_B = {
	0x4d940e93u, 0xddbcbd41u, 0x15ab8f92u, 0xf39789f5u,
	0xcf6509a7u, 0x4d5a9e4bu, 0x9d9f5e34u, 0x28e9fa9eu,
};


/**
 * @brief SM2 base point G in Jacobian coordinates
 */
const SM2_JACOBIAN_POINT SM2_GROUP = {
	{
		0x334c74c7u, 0x715a4589u, 0xf2660be1u, 0x8fe30bbfu,
		0x6a39c994u, 0x5f990446u, 0x1f198119u, 0x32c4ae2cu,
	},
	{
		0x2139f0a0u, 0x02df32e5u, 0xc62a4740u, 0xd0a9877cu,
		0x6b692153u, 0x59bdcee3u, 0xf4f6779cu, 0xbc3736a2u,
	},
	{
		1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u,
	},
};

/**
 * @brief Pointer to SM2 base point G
 */
const SM2_JACOBIAN_POINT* SM2_G = &SM2_GROUP;


/**
 * @brief SM2 curve order n
 */
const SM2_BN SM2_N = {
    0x39d54123u, 0x53bbf409u, 0x21c6052bu, 0x7203df6bu,
    0xffffffffu, 0xffffffffu, 0xffffffffu, 0xfffffffeu,
};


/**
 * @brief Constant u + 1 used in SM2 calculations
 */
const SM2_BN SM2_U_PLUS_ONE = {
    0x00000000u, 0x40000000u, 0xc0000000u, 0xffffffffu,
    0xffffffffu, 0xffffffffu, 0xbfffffffu, 0x3fffffffu,
};


/**
 * @brief Constant value 1 as SM2_BN
 */
const SM2_BN SM2_ONE = { 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u };


/**
 * @brief Constant value 2 as SM2_BN
 */
const SM2_BN SM2_TWO = { 2u, 0u, 0u, 0u, 0u, 0u, 0u, 0u };


/**
 * @brief Constant value 3 as SM2_BN
 */
const SM2_BN SM2_THREE = { 3u, 0u, 0u, 0u, 0u, 0u, 0u, 0u };


/**
 * @brief Check if SM2_BN value is valid (all words <= 0xffffffff)
 * @param a SM2_BN value to check
 * @return 1 if valid, -1 if invalid
 */
int XMen_lib_sm2_bn_check(const SM2_BN a)
{
    int err = 0;
    int i;
    for (i = 0; i < 8; i++) {
        if (a[i] > 0xffffffffu) {
            err++;
        }
    }
    if (err)
        return -1;
    else
        return 1;
}

/**
 * @brief Check if SM2_BN value is zero
 * @param a SM2_BN value to check
 * @return 1 if zero, 0 if non-zero
 */
int XMen_lib_sm2_bn_is_zero(const SM2_BN a)
{
    int i;
    for (i = 0; i < 8; i++) {
        if (a[i] != 0)
            return 0;
    }
    return 1;
}

/**
 * @brief Check if SM2_BN value is one
 * @param a SM2_BN value to check
 * @return 1 if one, 0 if not one
 */
int XMen_lib_sm2_bn_is_one(const SM2_BN a)
{
    int i;
    int ret = 1;

    for (i = 0; i < 7; i++) {
        if ((a[i + 1] != 0) || (a[0] != 1))
        {
            ret = 0;
            break;
        }

    }

    return ret;
}


/**
 * @brief Convert SM2_BN to byte array
 * @param a Input SM2_BN value
 * @param out Output buffer for 32-byte array
 */
void XMen_lib_sm2_bn_to_bytes(const SM2_BN a, uint8_t out[32])
{
    int i;
    uint32_t ai;
    for (i = 7; i >= 0; i--)
    {
        ai = (uint32_t)a[i];
        PUTU32(out, ai);
        out += SM2_uint32_SIZE;
    }
}


/**
 * @brief Convert byte array to SM2_BN
 * @param r Output SM2_BN value
 * @param in Input 32-byte array
 */
void XMen_lib_sm2_bn_from_bytes(SM2_BN r, const uint8_t in[32])
{
    int i;
    for (i = 7; i >= 0; i--) {
        r[i] = GETU32(in);
        in += SM2_uint32_SIZE;
    }
}


/**
 * @brief Convert ASN.1 integer to SM2_BN
 * @param r Output SM2_BN value
 * @param d Input ASN.1 integer buffer
 * @param dlen Length of input buffer
 * @return 1 on success, -1 on error
 */
int XMen_lib_sm2_bn_from_asn1_integer(SM2_BN r, const uint8_t* d, uint32_t dlen)
{
    uint8_t buf[32] = { 0 };
    if (!d || dlen == 0) {
        return -1;
    }
    if (dlen > sizeof(buf)) {
        return -1;
    }
    memcpy(buf + sizeof(buf) - dlen, d, dlen);
    XMen_lib_sm2_bn_from_bytes(r, buf);
    return 1;
}


/**
 * @brief Convert SM2_BN to binary string representation
 * @param a Input SM2_BN value to convert
 * @param bits Output buffer for binary string (256 chars)
 */
void XMen_lib_sm2_bn_to_bits(const SM2_BN a, char bits[256])
{
    int i, j;
    uint64_t w;
    // Process each 32-bit word from most significant to least
    for (i = 7; i >= 0; i--) {
        w = a[i];
        // Convert each bit to '0' or '1' character
        for (j = 0; j < 32; j++) {
            *bits++ = (w & 0x80000000) ? '1' : '0';
            w <<= 1;
        }
    }
}


/**
 * @brief Compare two SM2_BN values
 * @param a First SM2_BN value
 * @param b Second SM2_BN value
 * @return 1 if a > b, -1 if a < b, 0 if a == b
 */
int XMen_lib_sm2_bn_cmp(const SM2_BN a, const SM2_BN b)
{
    int i;
    // Compare words from most significant to least
    for (i = 7; i >= 0; i--) {
        if (a[i] > b[i])
            return 1;
        if (a[i] < b[i])
            return -1;
    }
    return 0;
}


/**
 * @brief Check if SM2_BN value is odd
 * @param a Input SM2_BN value to check
 * @return 1 if odd, 0 if even
 */
int XMen_lib_sm2_bn_is_odd(const SM2_BN a)
{
    return a[0] & 0x01;
}


/**
 * @brief Set SM2_BN to a 32-bit word value
 * @param r Output SM2_BN value
 * @param a 32-bit word value to set
 */
void XMen_lib_sm2_bn_set_word(SM2_BN r, uint32_t a)
{
    int i;
    // Set least significant word
    r[0] = a;
    // Clear all other words
    for (i = 1; i < 8; i++) {
        r[i] = 0;
    }
}


/**
 * @brief Right shift SM2_BN by specified number of bits
 * @param ret Output SM2_BN value
 * @param a Input SM2_BN value
 * @param nbits Number of bits to shift (0-31)
 * @return 1 on success, -1 if nbits > 31
 */
int XMen_lib_sm2_bn_rshift(SM2_BN ret, const SM2_BN a, unsigned int nbits)
{
    SM2_BN r;
    int i;

    if (nbits > 31) {
        return -1;
    }
    if (nbits == 0) {
        XMen_lib_sm2_bn_copy(ret, a);
    }

    // Perform right shift across words
    for (i = 0; i < 7; i++) {
        r[i] = a[i] >> nbits;
        r[i] |= (a[i + 1] << (32 - nbits)) & 0xffffffffu;
    }
    r[i] = a[i] >> nbits;
    XMen_lib_sm2_bn_copy(ret, r);
    return 1;
}


/**
 * @brief Add two SM2_BN values
 * @param r Output sum (r = a + b)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_bn_add(SM2_BN r, const SM2_BN a, const SM2_BN b)
{
    int i;
    // Add least significant words
    r[0] = a[0] + b[0];

    // Add remaining words with carry
    for (i = 1; i < 8; i++) {
        r[i] = a[i] + b[i] + (r[i - 1] >> 32);
    }
    // Mask to 32 bits
    for (i = 0; i < 7; i++) {
        r[i] &= 0xffffffffu;
    }
}


/**
 * @brief Subtract two SM2_BN values
 * @param ret Output difference (ret = a - b)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_bn_sub(SM2_BN ret, const SM2_BN a, const SM2_BN b)
{
    int i;
    SM2_BN r;
    // Subtract with borrow
    r[0] = ((uint64_t)1 << 32) + a[0] - b[0];
    for (i = 1; i < 7; i++) {
        r[i] = 0xffffffffu + a[i] - b[i] + (r[i - 1] >> 32);
        r[i - 1] &= 0xffffffffu;
    }
    r[i] = a[i] - b[i] + (r[i - 1] >> 32) - 1;
    r[i - 1] &= 0xffffffffu;
    XMen_lib_sm2_bn_copy(ret, r);
}


/**
 * @brief Generate random SM2_BN in range [0, range)
 * @param r Output random value
 * @param range Upper bound (exclusive)
 * @return 1 on success
 */
int XMen_lib_sm2_bn_rand_range(SM2_BN r, const SM2_BN range)
{
    uint8_t buf[32];
    // Generate random values until we get one in range
    do {
        XMen_lib_sm2_bn_from_bytes(r, buf);
    } while (XMen_lib_sm2_bn_cmp(r, range) >= 0);
    return 1;
}


/**
 * @brief Generate random SM2_BN modulo prime field p
 * @param r Output random value
 * @return 1 on success, -1 on error
 */
int XMen_lib_sm2_fp_rand(SM2_BN r)
{
    if (XMen_lib_sm2_bn_rand_range(r, SM2_P) != 1) {
        return -1;
    }
    return 1;
}


/**
 * @brief Add two SM2_BN values modulo prime field p
 * @param r Output sum (r = a + b mod p)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_fp_add(SM2_BN r, const SM2_BN a, const SM2_BN b)
{
    XMen_lib_sm2_bn_add(r, a, b);
    // Reduce modulo p if result >= p
    if (XMen_lib_sm2_bn_cmp(r, SM2_P) >= 0) {
        XMen_lib_sm2_bn_sub(r, r, SM2_P);
    }
}

/**
 * @brief Subtract two SM2_BN values modulo prime field p
 * @param r Output difference (r = a - b mod p)
 * @param a First operand (minuend)
 * @param b Second operand (subtrahend)
 */
void XMen_lib_sm2_fp_sub(SM2_BN r, const SM2_BN a, const SM2_BN b)
{
    // If a >= b, perform direct subtraction
    if (XMen_lib_sm2_bn_cmp(a, b) >= 0) {
        XMen_lib_sm2_bn_sub(r, a, b);
    }
    // If a < b, add p to a first then subtract b
    else {
        SM2_BN t;
        XMen_lib_sm2_bn_sub(t, SM2_P, b);  // t = p - b
        XMen_lib_sm2_bn_add(r, t, a);      // r = t + a = p - b + a
    }
}

/**
 * @brief Double a SM2_BN value modulo prime field p
 * @param r Output doubled value (r = 2a mod p)
 * @param a Input value to double
 */
void XMen_lib_sm2_fp_dbl(SM2_BN r, const SM2_BN a)
{
    XMen_lib_sm2_fp_add(r, a, a);  // r = a + a mod p
}

/**
 * @brief Triple a SM2_BN value modulo prime field p
 * @param r Output tripled value (r = 3a mod p)
 * @param a Input value to triple
 */
void XMen_lib_sm2_fp_tri(SM2_BN r, const SM2_BN a)
{
    SM2_BN t;
    XMen_lib_sm2_fp_dbl(t, a);     // t = 2a
    XMen_lib_sm2_fp_add(r, t, a);  // r = t + a = 3a mod p
}

/**
 * @brief Divide SM2_BN value by 2 modulo prime field p
 * @param r Output halved value (r = a/2 mod p)
 * @param a Input value to halve
 */
void XMen_lib_sm2_fp_div2(SM2_BN r, const SM2_BN a)
{
    int i;
    XMen_lib_sm2_bn_copy(r, a);
    // If input is odd, add p before dividing
    if (r[0] & 0x01) {
        XMen_lib_sm2_bn_add(r, r, SM2_P);
    }
    // Perform right shift by 1 bit
    for (i = 0; i < 7; i++) {
        r[i] = (r[i] >> 1) | ((r[i + 1] & 0x01) << 31);
    }
    r[i] >>= 1;
}

/**
 * @brief Negate SM2_BN value modulo prime field p
 * @param r Output negated value (r = -a mod p)
 * @param a Input value to negate
 */
void XMen_lib_sm2_fp_neg(SM2_BN r, const SM2_BN a)
{
    if (XMen_lib_sm2_bn_is_zero(a)) {
        XMen_lib_sm2_bn_copy(r, a);  // If a is 0, result is also 0
    }
    else {
        XMen_lib_sm2_bn_sub(r, SM2_P, a);  // r = p - a
    }
}

/**
 * @brief Multiply two SM2_BN values modulo prime field p
 * @param r Output product (r = a * b mod p)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_fp_mul(SM2_BN r, const SM2_BN a, const SM2_BN b)
{
    int i, j;
    uint64_t s[16] = { 0 };  // Intermediate results array
    SM2_BN d = { 0 };        // Temporary storage
    uint64_t u;

    // Compute raw product without modular reduction
    for (i = 0; i < 8; i++) {
        u = 0;
        for (j = 0; j < 8; j++) {
            u = s[i + j] + a[i] * b[j] + u;
            s[i + j] = u & 0xffffffffu;
            u >>= 32;
        }
        s[i + 8] = u;
    }

    // Fast modular reduction algorithm specific to SM2's prime p
    r[0] = s[0] + s[8] + s[9] + s[10] + s[11] + s[12] + ((s[13] + s[14] + s[15]) << 1);
    r[1] = s[1] + s[9] + s[10] + s[11] + s[12] + s[13] + ((s[14] + s[15]) << 1);
    r[2] = s[2];
    r[3] = s[3] + s[8] + s[11] + s[12] + s[14] + s[15] + (s[13] << 1);
    r[4] = s[4] + s[9] + s[12] + s[13] + s[15] + (s[14] << 1);
    r[5] = s[5] + s[10] + s[13] + s[14] + (s[15] << 1);
    r[6] = s[6] + s[11] + s[14] + s[15];
    r[7] = s[7] + s[8] + s[9] + s[10] + s[11] + s[15] + ((s[12] + s[13] + s[14] + s[15]) << 1);

    // Propagate carries
    for (i = 1; i < 8; i++) {
        r[i] += r[i - 1] >> 32;
        r[i - 1] &= 0xffffffffu;
    }

    // Final reduction step
    d[2] = s[8] + s[9] + s[13] + s[14];
    d[3] = d[2] >> 32;
    d[2] &= 0xffffffffu;
    XMen_lib_sm2_bn_sub(r, r, d);

    // Ensure result is less than p
    while (XMen_lib_sm2_bn_cmp(r, SM2_P) >= 0) {
        XMen_lib_sm2_bn_sub(r, r, SM2_P);
    }
}

/**
 * @brief Square a SM2_BN value modulo prime field p
 * @param r Output squared value (r = a^2 mod p)
 * @param a Input value to square
 */
void XMen_lib_sm2_fp_sqr(SM2_BN r, const SM2_BN a)
{
    XMen_lib_sm2_fp_mul(r, a, a);  // Use multiplication for squaring
}

/**
 * @brief Compute modular exponentiation in prime field p
 * @param r Output result (r = a^e mod p)
 * @param a Base value
 * @param e Exponent value
 */
void XMen_lib_sm2_fp_exp(SM2_BN r, const SM2_BN a, const SM2_BN e)
{
    SM2_BN t;
    uint32_t w;
    int i, j;

    XMen_lib_sm2_bn_set_one(t);  // Initialize result to 1

    // Square-and-multiply algorithm
    for (i = 7; i >= 0; i--) {
        w = (uint32_t)e[i];
        for (j = 0; j < 32; j++) {
            XMen_lib_sm2_fp_sqr(t, t);  // Square
            if (w & 0x80000000)             // If current bit is 1
                XMen_lib_sm2_fp_mul(t, t, a);  // Multiply by a
            w <<= 1;  // Move to next bit
        }
    }

    XMen_lib_sm2_bn_copy(r, t);  // Copy final result
}

/**
 * @brief Compute modular multiplicative inverse in prime field p
 * @param r Output inverse value (r = a^(-1) mod p)
 * @param a Input value to invert
 */
void XMen_lib_sm2_fp_inv(SM2_BN r, const SM2_BN a)
{
    SM2_BN a1 = { 0 };
    SM2_BN a2 = { 0 };
    SM2_BN a3 = { 0 };
    SM2_BN a4 = { 0 };
    SM2_BN a5 = { 0 };
    int i;

    // Addition chain for computing modular inverse
    // Based on Fermat's little theorem: a^(p-1) ≡ 1 (mod p)
    // Therefore a^(p-2) ≡ a^(-1) (mod p)
    XMen_lib_sm2_fp_sqr(a1, a);
    XMen_lib_sm2_fp_mul(a2, a1, a);
    XMen_lib_sm2_fp_sqr(a3, a2);
    XMen_lib_sm2_fp_sqr(a3, a3);
    XMen_lib_sm2_fp_mul(a3, a3, a2);
    XMen_lib_sm2_fp_sqr(a4, a3);
    XMen_lib_sm2_fp_sqr(a4, a4);
    XMen_lib_sm2_fp_sqr(a4, a4);
    XMen_lib_sm2_fp_sqr(a4, a4);
    XMen_lib_sm2_fp_mul(a4, a4, a3);
    XMen_lib_sm2_fp_sqr(a5, a4);
    for (i = 1; i < 8; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a5, a5, a4);
    for (i = 0; i < 8; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a5, a5, a4);
    for (i = 0; i < 4; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a5, a5, a3);
    XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a5, a5, a2);
    XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a5, a5, a);
    XMen_lib_sm2_fp_sqr(a4, a5);
    XMen_lib_sm2_fp_mul(a3, a4, a1);
    XMen_lib_sm2_fp_sqr(a5, a4);
    for (i = 1; i < 31; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a4, a5, a4);
    XMen_lib_sm2_fp_sqr(a4, a4);
    XMen_lib_sm2_fp_mul(a4, a4, a);
    XMen_lib_sm2_fp_mul(a3, a4, a2);
    for (i = 0; i < 33; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a2, a5, a3);
    XMen_lib_sm2_fp_mul(a3, a2, a3);
    for (i = 0; i < 32; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a2, a5, a3);
    XMen_lib_sm2_fp_mul(a3, a2, a3);
    XMen_lib_sm2_fp_mul(a4, a2, a4);
    for (i = 0; i < 32; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a2, a5, a3);
    XMen_lib_sm2_fp_mul(a3, a2, a3);
    XMen_lib_sm2_fp_mul(a4, a2, a4);
    for (i = 0; i < 32; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a2, a5, a3);
    XMen_lib_sm2_fp_mul(a3, a2, a3);
    XMen_lib_sm2_fp_mul(a4, a2, a4);
    for (i = 0; i < 32; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(a2, a5, a3);
    XMen_lib_sm2_fp_mul(a3, a2, a3);
    XMen_lib_sm2_fp_mul(a4, a2, a4);
    for (i = 0; i < 32; i++)
        XMen_lib_sm2_fp_sqr(a5, a5);
    XMen_lib_sm2_fp_mul(r, a4, a5);

    // Clean up temporary variables
    XMen_lib_sm2_bn_clean(a1);
    XMen_lib_sm2_bn_clean(a2);
    XMen_lib_sm2_bn_clean(a3);
    XMen_lib_sm2_bn_clean(a4);
    XMen_lib_sm2_bn_clean(a5);
}


/**
 * @brief Calculate square root in prime field p
 * @param r Output square root value (r = sqrt(a) mod p)
 * @param a Input value to calculate square root
 * @return 1 on success, -1 if a is not a quadratic residue
 */
int XMen_lib_sm2_fp_sqrt(SM2_BN r, const SM2_BN a)
{
    SM2_BN u;
    SM2_BN y;

    // Calculate exponent u = (p + 1)/4
    XMen_lib_sm2_bn_add(u, SM2_P, SM2_ONE);
    XMen_lib_sm2_bn_rshift(u, u, 2);

    // Calculate y = a^u mod p
    XMen_lib_sm2_fp_exp(y, a, u);

    // Verify y^2 = a mod p
    XMen_lib_sm2_fp_sqr(u, y);
    if (XMen_lib_sm2_bn_cmp(u, a) != 0) {
        return -1; // Not a quadratic residue
    }

    XMen_lib_sm2_bn_copy(r, y); // Copy result
    return 1;
}


/**
 * @brief Add two values modulo curve order n
 * @param r Output sum (r = a + b mod n)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_fn_add(SM2_Fn r, const SM2_Fn a, const SM2_Fn b)
{
    XMen_lib_sm2_bn_add(r, a, b);
    // Reduce modulo n if result >= n
    if (XMen_lib_sm2_bn_cmp(r, SM2_N) >= 0) {
        XMen_lib_sm2_bn_sub(r, r, SM2_N);
    }
}


/**
 * @brief Subtract two values modulo curve order n
 * @param r Output difference (r = a - b mod n)
 * @param a First operand (minuend)
 * @param b Second operand (subtrahend)
 */
void XMen_lib_sm2_fn_sub(SM2_Fn r, const SM2_Fn a, const SM2_Fn b)
{
    // If a >= b, perform direct subtraction
    if (XMen_lib_sm2_bn_cmp(a, b) >= 0) {
        XMen_lib_sm2_bn_sub(r, a, b);
    }
    // If a < b, add n to a first then subtract b
    else {
        SM2_BN t;
        XMen_lib_sm2_bn_add(t, a, SM2_N);
        XMen_lib_sm2_bn_sub(r, t, b);
    }
}


/**
 * @brief Negate value modulo curve order n
 * @param r Output negation (r = -a mod n)
 * @param a Input value to negate
 */
void XMen_lib_sm2_fn_neg(SM2_Fn r, const SM2_Fn a)
{
    if (XMen_lib_sm2_bn_is_zero(a)) {
        XMen_lib_sm2_bn_copy(r, a);
    }
    else {
        XMen_lib_sm2_bn_sub(r, SM2_N, a);
    }
}


/**
 * @brief Compare two 288-bit numbers
 * @param a First 288-bit number (9 words)
 * @param b Second 288-bit number (9 words)
 * @return 1 if a > b, -1 if a < b, 0 if a == b
 */
static int XMen_lib_sm2_bn288_cmp(const uint64_t a[9], const uint64_t b[9])
{
    int i;
    for (i = 8; i >= 0; i--) {
        if (a[i] > b[i])
            return 1;
        if (a[i] < b[i])
            return -1;
    }
    return 0;
}


/**
 * @brief Add two 288-bit numbers
 * @param r Output sum (9 words)
 * @param a First operand (9 words)
 * @param b Second operand (9 words)
 */
static void XMen_lib_sm2_bn288_add(uint64_t r[9], const uint64_t a[9], const uint64_t b[9])
{
    int i;
    // Add with carry propagation
    r[0] = a[0] + b[0];
    for (i = 1; i < 9; i++) {
        r[i] = a[i] + b[i] + (r[i - 1] >> 32);
    }
    // Mask to 32 bits per word
    for (i = 0; i < 8; i++) {
        r[i] &= 0xffffffffu;
    }
}


/**
 * @brief Subtract two 288-bit numbers
 * @param ret Output difference (9 words)
 * @param a First operand (minuend, 9 words)
 * @param b Second operand (subtrahend, 9 words)
 */
static void XMen_lib_sm2_bn288_sub(uint64_t ret[9], const uint64_t a[9], const uint64_t b[9])
{
    int i;
    uint64_t r[9];

    // Subtract with borrow propagation
    r[0] = ((uint64_t)1 << 32) + a[0] - b[0];
    for (i = 1; i < 8; i++) {
        r[i] = 0xffffffffu + a[i] - b[i] + (r[i - 1] >> 32);
        r[i - 1] &= 0xffffffffu;
    }
    r[i] = a[i] - b[i] + (r[i - 1] >> 32) - 1;
    r[i - 1] &= 0xffffffffu;

    // Copy result
    for (i = 0; i < 9; i++) {
        ret[i] = r[i];
    }
}

/**
 * @brief Multiply two values modulo curve order n using Barrett reduction
 * @param ret Output product (ret = a * b mod n)
 * @param a First operand
 * @param b Second operand
 */
void XMen_lib_sm2_fn_mul(SM2_BN ret, const SM2_BN a, const SM2_BN b)
{
    SM2_BN r;
    // Barrett reduction constant mu
    static const uint64_t mu[9] = {
        0xf15149a0u, 0x12ac6361u, 0xfa323c01u, 0x8dfc2096u, 1u, 1u, 1u, 1u, 1u,
    };

    uint64_t s[18];
    uint64_t zh[9];
    uint64_t zl[9];
    uint64_t q[9];
    uint64_t w;
    int i, j;

    // Compute product s = a * b
    for (i = 0; i < 8; i++) {
        s[i] = 0;
    }
    for (i = 0; i < 8; i++) {
        w = 0;
        for (j = 0; j < 8; j++) {
            w += s[i + j] + a[i] * b[j];
            s[i + j] = w & 0xffffffffu;
            w >>= 32;
        }
        s[i + 8] = w;
    }

    // Split s into high and low parts
    for (i = 0; i < 9; i++) {
        zl[i] = s[i];
        zh[i] = s[7 + i];
    }

    // Compute quotient estimate q = floor(zh * mu / 2^288)
    for (i = 0; i < 9; i++) {
        s[i] = 0;
    }
    for (i = 0; i < 9; i++) {
        w = 0;
        for (j = 0; j < 9; j++) {
            w += s[i + j] + zh[i] * mu[j];
            s[i + j] = w & 0xffffffffu;
            w >>= 32;
        }
        s[i + 9] = w;
    }
    for (i = 0; i < 8; i++) {
        q[i] = s[9 + i];
    }

    // Compute q * n
    for (i = 0; i < 17; i++) {
        s[i] = 0;
    }
    for (i = 0; i < 8; i++) {
        w = 0;
        for (j = 0; j < 8; j++) {
            w += s[i + j] + q[i] * SM2_N[j];
            s[i + j] = w & 0xffffffffu;
            w >>= 32;
        }
        s[i + 8] = w;
    }
    for (i = 0; i < 9; i++) {
        q[i] = s[i];
    }

    // Compute final result r = zl - q * n
    if (XMen_lib_sm2_bn288_cmp(zl, q)) {
        XMen_lib_sm2_bn288_sub(zl, zl, q);
    }
    else {
        uint64_t c[9] = { 0,0,0,0,0,0,0,0,0x100000000 };
        XMen_lib_sm2_bn288_sub(q, c, q);
        XMen_lib_sm2_bn288_add(zl, q, zl);
    }
    for (i = 0; i < 8; i++) {
        r[i] = zl[i];
    }
    r[7] += zl[8] << 32;

    // Final reduction modulo n
    while (XMen_lib_sm2_bn_cmp(r, SM2_N) >= 0) {
        XMen_lib_sm2_bn_sub(r, r, SM2_N);
    }
    XMen_lib_sm2_bn_copy(ret, r);
}

/**
 * @brief Multiply a field element by a word modulo curve order n
 * @param r Output product (r = a * b mod n)
 * @param a Input field element
 * @param b Input word to multiply by
 */
void XMen_lib_sm2_fn_mul_word(SM2_Fn r, const SM2_Fn a, uint32_t b)
{
    SM2_Fn t;
    XMen_lib_sm2_bn_set_word(t, b);  // Convert word to field element
    XMen_lib_sm2_fn_mul(r, a, t);    // Multiply field elements
}


/**
 * @brief Square a field element modulo curve order n
 * @param r Output square (r = a^2 mod n)
 * @param a Input value to square
 */
void XMen_lib_sm2_fn_sqr(SM2_BN r, const SM2_BN a)
{
    XMen_lib_sm2_fn_mul(r, a, a);  // Square is multiplication by itself
}


/**
 * @brief Modular exponentiation modulo curve order n
 * @param r Output result (r = a^e mod n)
 * @param a Base value
 * @param e Exponent
 */
void XMen_lib_sm2_fn_exp(SM2_BN r, const SM2_BN a, const SM2_BN e)
{
    SM2_BN t;
    uint32_t w;
    int i, j;

    // Initialize result to 1
    XMen_lib_sm2_bn_set_one(t);

    // Square-and-multiply algorithm
    for (i = 7; i >= 0; i--) {
        w = (uint32_t)e[i];
        for (j = 0; j < 32; j++) {
            XMen_lib_sm2_fn_sqr(t, t);  // Square
            if (w & 0x80000000) {           // Multiply if bit is 1
                XMen_lib_sm2_fn_mul(t, t, a);
            }
            w <<= 1;
        }
    }
    XMen_lib_sm2_bn_copy(r, t);
}


/**
 * @brief Compute modular multiplicative inverse modulo curve order n
 * @param r Output inverse (r = a^(-1) mod n)
 * @param a Input value to invert
 */
void XMen_lib_sm2_fn_inv(SM2_BN r, const SM2_BN a)
{
    SM2_BN e;
    // Use Fermat's little theorem: a^(-1) = a^(n-2) mod n
    XMen_lib_sm2_bn_sub(e, SM2_N, SM2_TWO);
    XMen_lib_sm2_fn_exp(r, a, e);
}


/**
 * @brief Generate random field element modulo curve order n
 * @param r Output random value
 * @return 1 on success, -1 on error
 */
int XMen_lib_sm2_fn_rand(SM2_BN r)
{
    if (XMen_lib_sm2_bn_rand_range(r, SM2_N) != 1) {
        return -1;
    }
    return 1;
}


/**
 * @brief Initialize a Jacobian point to (1:1:0)
 * @param R Point to initialize
 */
void XMen_lib_sm2_jacobian_point_init(SM2_JACOBIAN_POINT* R)
{
    memset(R, 0, sizeof(SM2_JACOBIAN_POINT));
    R->X[0] = 1;
    R->Y[0] = 1;
}


/**
 * @brief Check if a Jacobian point is the point at infinity
 * @param P Point to check
 * @return 1 if point is at infinity, 0 otherwise
 */
int XMen_lib_sm2_jacobian_point_is_at_infinity(const SM2_JACOBIAN_POINT* P)
{
    return XMen_lib_sm2_bn_is_zero(P->Z);
}


/**
 * @brief Set Jacobian point coordinates from affine coordinates
 * @param R Output point in Jacobian coordinates
 * @param x Affine x-coordinate
 * @param y Affine y-coordinate
 */
void XMen_lib_sm2_jacobian_point_set_xy(SM2_JACOBIAN_POINT* R, const SM2_BN x, const SM2_BN y)
{
    XMen_lib_sm2_bn_copy(R->X, x);
    XMen_lib_sm2_bn_copy(R->Y, y);
    XMen_lib_sm2_bn_set_one(R->Z);  // Set Z coordinate to 1
}


/**
 * @brief Get affine coordinates from Jacobian point
 * @param P Input point in Jacobian coordinates
 * @param x Output affine x-coordinate
 * @param y Output affine y-coordinate (can be NULL if not needed)
 */
void XMen_lib_sm2_jacobian_point_get_xy(const SM2_JACOBIAN_POINT* P, SM2_BN x, SM2_BN y)
{
    if (XMen_lib_sm2_bn_is_one(P->Z)) {
        // If Z=1, coordinates are already in affine form
        XMen_lib_sm2_bn_copy(x, P->X);
        if (y) {
            XMen_lib_sm2_bn_copy(y, P->Y);
        }
    }
    else {
        // Convert from Jacobian to affine coordinates
        SM2_BN z_inv;
        XMen_lib_sm2_fp_inv(z_inv, P->Z);
        if (y) {
            XMen_lib_sm2_fp_mul(y, P->Y, z_inv);
        }
        XMen_lib_sm2_fp_sqr(z_inv, z_inv);
        XMen_lib_sm2_fp_mul(x, P->X, z_inv);
        if (y) {
            XMen_lib_sm2_fp_mul(y, y, z_inv);
        }
    }
}


/**
 * @brief Check if a point lies on the SM2 curve
 * @param P Point to check
 * @return 1 if point is on curve, -1 if not
 */
int XMen_lib_sm2_jacobian_point_is_on_curve(const SM2_JACOBIAN_POINT* P)
{
    SM2_BN t0, t1, t2;

    // Check curve equation y^2 = x^3 + ax + b
    if (XMen_lib_sm2_bn_is_one(P->Z)) {
        // For Z=1, check y^2 = x^3 + ax + b directly
        XMen_lib_sm2_fp_sqr(t0, P->Y);
        XMen_lib_sm2_fp_add(t0, t0, P->X);
        XMen_lib_sm2_fp_add(t0, t0, P->X);
        XMen_lib_sm2_fp_add(t0, t0, P->X);
        XMen_lib_sm2_fp_sqr(t1, P->X);
        XMen_lib_sm2_fp_mul(t1, t1, P->X);
        XMen_lib_sm2_fp_add(t1, t1, SM2_B);
    }
    else {
        // For Z≠1, check in Jacobian coordinates
        XMen_lib_sm2_fp_sqr(t0, P->Y);
        XMen_lib_sm2_fp_sqr(t1, P->Z);
        XMen_lib_sm2_fp_sqr(t2, t1);
        XMen_lib_sm2_fp_mul(t1, t1, t2);
        XMen_lib_sm2_fp_mul(t1, t1, SM2_B);
        XMen_lib_sm2_fp_mul(t2, t2, P->X);
        XMen_lib_sm2_fp_add(t0, t0, t2);
        XMen_lib_sm2_fp_add(t0, t0, t2);
        XMen_lib_sm2_fp_add(t0, t0, t2);
        XMen_lib_sm2_fp_sqr(t2, P->X);
        XMen_lib_sm2_fp_mul(t2, t2, P->X);
        XMen_lib_sm2_fp_add(t1, t1, t2);
    }

    return (XMen_lib_sm2_bn_cmp(t0, t1) == 0) ? 1 : -1;
}


/**
 * @brief Negate a point in Jacobian coordinates
 * @param R Output negated point
 * @param P Input point to negate
 */
void XMen_lib_sm2_jacobian_point_neg(SM2_JACOBIAN_POINT* R, const SM2_JACOBIAN_POINT* P)
{
    XMen_lib_sm2_bn_copy(R->X, P->X);
    XMen_lib_sm2_fp_neg(R->Y, P->Y);  // Negate Y coordinate
    XMen_lib_sm2_bn_copy(R->Z, P->Z);
}


/**
 * @brief Double a point in Jacobian coordinates
 * @param R Output doubled point
 * @param P Input point to double
 */
void XMen_lib_sm2_jacobian_point_dbl(SM2_JACOBIAN_POINT* R, const SM2_JACOBIAN_POINT* P)
{
    const uint64_t* X1 = P->X;
    const uint64_t* Y1 = P->Y;
    const uint64_t* Z1 = P->Z;
    SM2_BN T1, T2, T3;
    SM2_BN X3, Y3, Z3;

    // Handle point at infinity
    if (XMen_lib_sm2_jacobian_point_is_at_infinity(P)) {
        XMen_lib_sm2_jacobian_point_copy(R, P);
        return;
    }

    // Point doubling formula for Jacobian coordinates
    XMen_lib_sm2_fp_sqr(T1, Z1);          // T1 = Z1^2
    XMen_lib_sm2_fp_sub(T2, X1, T1);      // T2 = X1 - T1
    XMen_lib_sm2_fp_add(T1, X1, T1);      // T1 = X1 + T1
    XMen_lib_sm2_fp_mul(T2, T2, T1);      // T2 = T2 * T1
    XMen_lib_sm2_fp_tri(T2, T2);          // T2 = 3 * T2
    XMen_lib_sm2_fp_dbl(Y3, Y1);          // Y3 = 2 * Y1
    XMen_lib_sm2_fp_mul(Z3, Y3, Z1);      // Z3 = Y3 * Z1
    XMen_lib_sm2_fp_sqr(Y3, Y3);          // Y3 = Y3^2
    XMen_lib_sm2_fp_mul(T3, Y3, X1);      // T3 = Y3 * X1
    XMen_lib_sm2_fp_sqr(Y3, Y3);          // Y3 = Y3^2
    XMen_lib_sm2_fp_div2(Y3, Y3);         // Y3 = Y3/2
    XMen_lib_sm2_fp_sqr(X3, T2);          // X3 = T2^2
    XMen_lib_sm2_fp_dbl(T1, T3);          // T1 = 2 * T3
    XMen_lib_sm2_fp_sub(X3, X3, T1);      // X3 = X3 - T1
    XMen_lib_sm2_fp_sub(T1, T3, X3);      // T1 = T3 - X3
    XMen_lib_sm2_fp_mul(T1, T1, T2);      // T1 = T1 * T2
    XMen_lib_sm2_fp_sub(Y3, T1, Y3);      // Y3 = T1 - Y3

    // Set result coordinates
    XMen_lib_sm2_bn_copy(R->X, X3);
    XMen_lib_sm2_bn_copy(R->Y, Y3);
    XMen_lib_sm2_bn_copy(R->Z, Z3);
}


/**
 * @brief Add two points in Jacobian coordinates
 * @param R Output point for sum (R = P + Q)
 * @param P First input point
 * @param Q Second input point (must be in affine coordinates with Z=1)
 */
void XMen_lib_sm2_jacobian_point_add(SM2_JACOBIAN_POINT* R, const SM2_JACOBIAN_POINT* P, const SM2_JACOBIAN_POINT* Q)
{
    const uint64_t* X1 = P->X;
    const uint64_t* Y1 = P->Y;
    const uint64_t* Z1 = P->Z;
    const uint64_t* x2 = Q->X;  // Q is in affine coordinates
    const uint64_t* y2 = Q->Y;
    SM2_BN T1;
    SM2_BN T2;
    SM2_BN T3;
    SM2_BN T4;
    SM2_BN X3;
    SM2_BN Y3;
    SM2_BN Z3;

    // Handle special cases - if either point is infinity
    if (XMen_lib_sm2_jacobian_point_is_at_infinity(Q)) {
        XMen_lib_sm2_jacobian_point_copy(R, P);
        return;
    }

    if (XMen_lib_sm2_jacobian_point_is_at_infinity(P)) {
        XMen_lib_sm2_jacobian_point_copy(R, Q);
        return;
    }

    // Q must be in affine coordinates (Z=1)
    assert(XMen_lib_sm2_bn_is_one(Q->Z));

    // Point addition formula for mixed Jacobian-affine coordinates
    XMen_lib_sm2_fp_sqr(T1, Z1);          // T1 = Z1^2
    XMen_lib_sm2_fp_mul(T2, T1, Z1);      // T2 = Z1^3
    XMen_lib_sm2_fp_mul(T1, T1, x2);      // T1 = X2*Z1^2
    XMen_lib_sm2_fp_mul(T2, T2, y2);      // T2 = Y2*Z1^3
    XMen_lib_sm2_fp_sub(T1, T1, X1);      // T1 = X2*Z1^2 - X1
    XMen_lib_sm2_fp_sub(T2, T2, Y1);      // T2 = Y2*Z1^3 - Y1

    // Check if points are equal or opposite
    if (XMen_lib_sm2_bn_is_zero(T1)) {
        if (XMen_lib_sm2_bn_is_zero(T2)) {
            // Points are equal - use point doubling
            SM2_JACOBIAN_POINT _Q, * Q = &_Q;
            XMen_lib_sm2_jacobian_point_set_xy(Q, x2, y2);
            XMen_lib_sm2_jacobian_point_dbl(R, Q);
            return;
        }
        else {
            // Points are opposite - result is infinity
            XMen_lib_sm2_jacobian_point_set_infinity(R);
            return;
        }
    }

    // Complete point addition
    XMen_lib_sm2_fp_mul(Z3, Z1, T1);      // Z3 = Z1*T1
    XMen_lib_sm2_fp_sqr(T3, T1);          // T3 = T1^2
    XMen_lib_sm2_fp_mul(T4, T3, T1);      // T4 = T1^3
    XMen_lib_sm2_fp_mul(T3, T3, X1);      // T3 = X1*T1^2
    XMen_lib_sm2_fp_dbl(T1, T3);          // T1 = 2*X1*T1^2
    XMen_lib_sm2_fp_sqr(X3, T2);          // X3 = T2^2
    XMen_lib_sm2_fp_sub(X3, X3, T1);      // X3 = T2^2 - 2*X1*T1^2
    XMen_lib_sm2_fp_sub(X3, X3, T4);      // X3 = T2^2 - 2*X1*T1^2 - T1^3
    XMen_lib_sm2_fp_sub(T3, T3, X3);      // T3 = X1*T1^2 - X3
    XMen_lib_sm2_fp_mul(T3, T3, T2);      // T3 = T2*(X1*T1^2 - X3)
    XMen_lib_sm2_fp_mul(T4, T4, Y1);      // T4 = Y1*T1^3
    XMen_lib_sm2_fp_sub(Y3, T3, T4);      // Y3 = T2*(X1*T1^2 - X3) - Y1*T1^3

    // Set result coordinates
    XMen_lib_sm2_bn_copy(R->X, X3);
    XMen_lib_sm2_bn_copy(R->Y, Y3);
    XMen_lib_sm2_bn_copy(R->Z, Z3);
}


/**
 * @brief Subtract two points in Jacobian coordinates
 * @param R Output point for difference (R = P - Q)
 * @param P First input point
 * @param Q Second input point
 */
void XMen_lib_sm2_jacobian_point_sub(SM2_JACOBIAN_POINT* R, const SM2_JACOBIAN_POINT* P, const SM2_JACOBIAN_POINT* Q)
{
    SM2_JACOBIAN_POINT _T, * T = &_T;
    XMen_lib_sm2_jacobian_point_neg(T, Q);  // T = -Q
    XMen_lib_sm2_jacobian_point_add(R, P, T);  // R = P + (-Q)
}


/**
 * @brief Multiply point by scalar in Jacobian coordinates
 * @param R Output point for product (R = k*P)
 * @param k Scalar multiplier
 * @param P Input point to multiply
 */
void XMen_lib_sm2_jacobian_point_mul(SM2_JACOBIAN_POINT* R, const SM2_BN k, const SM2_JACOBIAN_POINT* P)
{
    char bits[257] = { 0 };
    SM2_JACOBIAN_POINT _Q, * Q = &_Q;
    SM2_JACOBIAN_POINT _T, * T = &_T;
    int i;

    // Convert P to affine coordinates if needed
    if (!XMen_lib_sm2_bn_is_one(P->Z)) {
        SM2_BN x;
        SM2_BN y;
        XMen_lib_sm2_jacobian_point_get_xy(P, x, y);
        XMen_lib_sm2_jacobian_point_set_xy(T, x, y);
        P = T;
    }

    // Double-and-add algorithm for point multiplication
    XMen_lib_sm2_jacobian_point_set_infinity(Q);  // Q = infinity
    XMen_lib_sm2_bn_to_bits(k, bits);  // Convert k to binary
    for (i = 0; i < 256; i++) {
        XMen_lib_sm2_jacobian_point_dbl(Q, Q);  // Q = 2Q
        if (bits[i] == '1') {
            XMen_lib_sm2_jacobian_point_add(Q, Q, P);  // Q = Q + P
        }
    }
    XMen_lib_sm2_jacobian_point_copy(R, Q);
}


/**
 * @brief Convert Jacobian point to compressed byte format
 * @param P Input point in Jacobian coordinates
 * @param out Output buffer for 64-byte point representation
 */
void XMen_lib_sm2_jacobian_point_to_bytes(const SM2_JACOBIAN_POINT* P, uint8_t out[64])
{
    SM2_BN x;
    SM2_BN y;
    XMen_lib_sm2_jacobian_point_get_xy(P, x, y);  // Convert to affine coordinates
    XMen_lib_sm2_bn_to_bytes(x, out);  // First 32 bytes = x coordinate
    XMen_lib_sm2_bn_to_bytes(y, out + 32);  // Second 32 bytes = y coordinate
}


/**
 * @brief Convert byte format to Jacobian point
 * @param P Output point in Jacobian coordinates
 * @param in Input 64-byte point representation
 */
void XMen_lib_sm2_jacobian_point_from_bytes(SM2_JACOBIAN_POINT* P, const uint8_t in[64])
{
    XMen_lib_sm2_bn_from_bytes(P->X, in);  // First 32 bytes = x coordinate
    XMen_lib_sm2_bn_from_bytes(P->Y, in + 32);  // Second 32 bytes = y coordinate
    XMen_lib_sm2_bn_set_word(P->Z, 1);  // Set Z=1 for affine coordinates
}


/**
 * @brief Multiply generator point by scalar
 * @param R Output point for product (R = k*G)
 * @param k Scalar multiplier
 */
void XMen_lib_sm2_jacobian_point_mul_generator(SM2_JACOBIAN_POINT* R, const SM2_BN k)
{
    XMen_lib_sm2_jacobian_point_mul(R, k, SM2_G);  // R = k*G
}


/**
 * @brief Compute linear combination of points (R = t*P + s*G)
 * @param R Output point for sum
 * @param t First scalar multiplier
 * @param P First input point
 * @param s Second scalar multiplier
 */
void XMen_lib_sm2_jacobian_point_mul_sum(SM2_JACOBIAN_POINT* R, const SM2_BN t, const SM2_JACOBIAN_POINT* P, const SM2_BN s)
{
    SM2_JACOBIAN_POINT _sG, * sG = &_sG;
    SM2_BN x;
    SM2_BN y;

    // Compute s*G
    XMen_lib_sm2_jacobian_point_mul_generator(sG, s);

    // Compute t*P and convert to affine coordinates
    XMen_lib_sm2_jacobian_point_mul(R, t, P);
    XMen_lib_sm2_jacobian_point_get_xy(R, x, y);
    XMen_lib_sm2_jacobian_point_set_xy(R, x, y);

    // Add results: R = s*G + t*P
    XMen_lib_sm2_jacobian_point_add(R, sG, R);
}

/**
 * @brief Check if point is on the SM2 curve
 * @param P Input point to check
 * @return 1 if point is on curve, 0 otherwise
 */
int XMen_lib_sm2_point_is_on_curve(const SM2_POINT* P)
{
    SM2_JACOBIAN_POINT T;
    XMen_lib_sm2_jacobian_point_from_bytes(&T, (const uint8_t*)P);
    return XMen_lib_sm2_jacobian_point_is_on_curve(&T);
}


/**
 * @brief Compute point from x-coordinate and y-bit
 * @param P Output point
 * @param x Input x-coordinate (32 bytes)
 * @param y Y-bit indicating which y value to use (0x02 or 0x03)
 * @return 1 on success, -1 on error
 */
int XMen_lib_sm2_point_from_x(SM2_POINT* P, const uint8_t x[32], int y)
{
    SM2_BN _x, _y, _g, _z;
    XMen_lib_sm2_bn_from_bytes(_x, x);

    // Compute right side of curve equation: g = x^3 - 3x + b
    XMen_lib_sm2_fp_sqr(_g, _x);
    XMen_lib_sm2_fp_sub(_g, _g, SM2_THREE);
    XMen_lib_sm2_fp_mul(_g, _g, _x);
    XMen_lib_sm2_fp_add(_g, _g, SM2_B);

    // Compute square root
    XMen_lib_sm2_fp_exp(_y, _g, SM2_U_PLUS_ONE);

    // Verify solution
    XMen_lib_sm2_fp_sqr(_z, _y);
    if (XMen_lib_sm2_bn_cmp(_z, _g)) {
        return -1;
    }

    // Select correct y value based on y-bit
    if ((y == 0x02 && XMen_lib_sm2_bn_is_odd(_y)) || ((y == 0x03) && !XMen_lib_sm2_bn_is_odd(_y))) {
        XMen_lib_sm2_fp_neg(_y, _y);
    }

    // Convert to bytes and clean up
    XMen_lib_sm2_bn_to_bytes(_x, P->x);
    XMen_lib_sm2_bn_to_bytes(_y, P->y);

    XMen_lib_sm2_bn_clean(_x);
    XMen_lib_sm2_bn_clean(_y);
    XMen_lib_sm2_bn_clean(_g);
    XMen_lib_sm2_bn_clean(_z);

    // Verify point is on curve
    if (XMen_lib_sm2_point_is_on_curve(P) != 1) {
        return -1;
    }
    return 1;
}


/**
 * @brief Create point from x and y coordinates
 * @param P Output point
 * @param x Input x-coordinate (32 bytes)
 * @param y Input y-coordinate (32 bytes)
 * @return 1 if point is on curve, -1 otherwise
 */
int XMen_lib_sm2_point_from_xy(SM2_POINT* P, const uint8_t x[32], const uint8_t y[32])
{
    memcpy(P->x, x, 32);
    memcpy(P->y, y, 32);
    return XMen_lib_sm2_point_is_on_curve(P);
}


/**
 * @brief Add two points in affine coordinates
 * @param R Output point for sum (R = P + Q)
 * @param P First input point
 * @param Q Second input point
 * @return 1 on success
 */
int XMen_lib_sm2_point_add(SM2_POINT* R, const SM2_POINT* P, const SM2_POINT* Q)
{
    SM2_JACOBIAN_POINT P_;
    SM2_JACOBIAN_POINT Q_;

    // Convert to Jacobian coordinates
    XMen_lib_sm2_jacobian_point_from_bytes(&P_, (uint8_t*)P);
    XMen_lib_sm2_jacobian_point_from_bytes(&Q_, (uint8_t*)Q);
    
    // Add points and convert back
    XMen_lib_sm2_jacobian_point_add(&P_, &P_, &Q_);
    XMen_lib_sm2_jacobian_point_to_bytes(&P_, (uint8_t*)R);

    return 1;
}


/**
 * @brief Subtract two points in affine coordinates
 * @param R Output point for difference (R = P - Q)
 * @param P First input point
 * @param Q Second input point
 * @return 1 on success
 */
int XMen_lib_sm2_point_sub(SM2_POINT* R, const SM2_POINT* P, const SM2_POINT* Q)
{
    SM2_JACOBIAN_POINT P_;
    SM2_JACOBIAN_POINT Q_;

    // Convert to Jacobian coordinates
    XMen_lib_sm2_jacobian_point_from_bytes(&P_, (uint8_t*)P);
    XMen_lib_sm2_jacobian_point_from_bytes(&Q_, (uint8_t*)Q);
    
    // Subtract points and convert back
    XMen_lib_sm2_jacobian_point_sub(&P_, &P_, &Q_);
    XMen_lib_sm2_jacobian_point_to_bytes(&P_, (uint8_t*)R);

    return 1;
}


/**
 * @brief Negate a point in affine coordinates
 * @param R Output point for negation (R = -P)
 * @param P Input point to negate
 * @return 1 on success
 */
int XMen_lib_sm2_point_neg(SM2_POINT* R, const SM2_POINT* P)
{
    SM2_JACOBIAN_POINT P_;

    // Convert to Jacobian coordinates
    XMen_lib_sm2_jacobian_point_from_bytes(&P_, (uint8_t*)P);
    
    // Negate point and convert back
    XMen_lib_sm2_jacobian_point_neg(&P_, &P_);
    XMen_lib_sm2_jacobian_point_to_bytes(&P_, (uint8_t*)R);

    return 1;
}


/**
 * @brief Double a point in affine coordinates
 * @param R Output point for doubling (R = 2P)
 * @param P Input point to double
 * @return 1 on success
 */
int XMen_lib_sm2_point_dbl(SM2_POINT* R, const SM2_POINT* P)
{
    SM2_JACOBIAN_POINT P_;

    // Convert to Jacobian coordinates
    XMen_lib_sm2_jacobian_point_from_bytes(&P_, (uint8_t*)P);
    
    // Double point and convert back
    XMen_lib_sm2_jacobian_point_dbl(&P_, &P_);
    XMen_lib_sm2_jacobian_point_to_bytes(&P_, (uint8_t*)R);

    return 1;
}


/**
 * @brief Multiply a point by a scalar in affine coordinates
 * @param R Output point for scalar multiplication (R = kP)
 * @param k Input scalar value as 32-byte array
 * @param P Input point to multiply
 * @return 1 on success
 */
int XMen_lib_sm2_point_mul(SM2_POINT* R, const uint8_t k[32], const SM2_POINT* P)
{
    SM2_BN _k;
    SM2_JACOBIAN_POINT _P;

    // Convert scalar from bytes to big number
    XMen_lib_sm2_bn_from_bytes(_k, k);
    
    // Convert point to Jacobian coordinates
    XMen_lib_sm2_jacobian_point_from_bytes(&_P, (uint8_t*)P);
    
    // Perform scalar multiplication
    XMen_lib_sm2_jacobian_point_mul(&_P, _k, &_P);
    
    // Convert result back to affine coordinates
    XMen_lib_sm2_jacobian_point_to_bytes(&_P, (uint8_t*)R);

    // Clean up temporary variables
    XMen_lib_sm2_bn_clean(_k);
    return 1;
}


/**
 * @brief Convert point from octet string to internal format
 * @param P Output point in internal format
 * @param in Input octet string buffer
 * @param inlen Length of input buffer
 * @return 1 on success, -1 on error
 *
 * Supports compressed (0x02/0x03 + 32 bytes) and 
 * uncompressed (0x04 + 64 bytes) point formats
 */
int XMen_lib_sm2_point_from_octets(SM2_POINT* P, const uint8_t* in, uint32_t inlen)
{
    // Handle compressed point format (33 bytes total)
    if ((*in == 0x02 || *in == 0x03) && inlen == 33) {
        // Reconstruct point from x-coordinate and compression byte
        if (XMen_lib_sm2_point_from_x(P, in + 1, *in) != 1) {
            return -1;
        }
    }
    // Handle uncompressed point format (65 bytes total) 
    else if (*in == 0x04 && inlen == 65) {
        // Parse x and y coordinates directly
        if (XMen_lib_sm2_point_from_xy(P, in + 1, in + 33) != 1) {
            return -1;
        }
    }
    else {
        return -1;  // Invalid format
    }
    return 1;
}



